========================================================================= 
Log Path: ./log/log_A5_2.log 
Program Path: /Users/james/Dropbox/Research/Diplomacy/Replication/A5_part_2.R 
Working Directory: /Users/james/Dropbox/Research/Diplomacy/Replication 
User Name: james 
R Version: 4.1.2 (2021-11-01) 
Machine: Jamess-MacBook-Air.local x86_64 
[1] "Operating System: Darwin 20.6.0 Darwin Kernel Version 20.6.0: Tue Feb 22 21:10:41 PST 2022; root:xnu-7195.141.26~1/RELEASE_X86_64"
[1] "Base Packages: stats graphics grDevices utils datasets methods base\nOther Packages: logr_1.3.0 gridExtra_2.3 ggpubr_0.4.0 estimatr_0.30.2 coefplot_1.2.7\n                rio_0.5.29 xtable_1.8-4 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7\n                purrr_0.3.4 readr_1.4.0 tidyr_1.1.3 tibble_3.1.2 ggplot2_3.3.6\n                tidyverse_1.3.1 "
Log Start Time: 2022-06-18 13:21:35 
========================================================================= 

> library(logr)
> log_open("log_A5_2.log")
> log_code()
> ####################
> #REPLICATION FILES: APPENDIX 5
> #Article: "When Does Online Public Diplomacy Succeed? Evidence from China's ‘Wolf Warrior’ Diplomats"
> #Authors: Daniel Mattingly and James Sundquist
> #This Version: June 18, 2022
> 
> 
> ####################
> #Description:
> # This script conducts our primary analysis (estimating Average Treatment Effects using
> # Ordinary Least Squares) on different versions of composite outcomes.
> # The datasets are included in the replication folder, but can also be reconstructed in the file
> # A5_part_1.R
> 
> #The script was written with R version 4.1.2
> #It was created and tested on Mac OS X (11.6.5).
> 
> 
> library(tidyverse)
> library(estimatr)
> library(coefplot)
> 
> setwd("~/Dropbox/Research/Diplomacy/Replication")
> 
> load("wave1_both_loadings.Rdata")
> load("wave2_both_loadings.Rdata")
> 
> may <- pre
> june <- post
> pooled <- rbind(may, june)
> 
> ############################################
> #Fit Models
> ############################################
> 
> # PCA fit on pooled data
> m.gov.may <- lm(china_gov_pca ~ t_, data = may)
> m.gov.june <- lm(china_gov_pca ~ t_, data = june)
> # PCA fit separately for each wave
> m.gov.may.waves <- lm(china_gov_pca_waves ~ t_, data = may)
> m.gov.june.waves <- lm(china_gov_pca_waves ~ t_, data = june)
> 
> m.people.may <- lm(china_people_pca ~ t_, data = may)
> m.people.june <- lm(china_people_pca ~ t_, data = june)
> m.people.may.waves <- lm(china_people_pca_waves ~ t_, data = may)
> m.people.june.waves <- lm(china_people_pca_waves ~ t_, data = june)
> 
> m.policy.may <- lm(china_policy_pca ~ t_, data = may)
> m.policy.june <- lm(china_policy_pca ~ t_, data = june)
> m.policy.may.waves <- lm(china_policy_pca_waves ~ t_, data = may)
> m.policy.june.waves <- lm(china_policy_pca_waves ~ t_, data = june)
> 
> m.covid.may <- lm(china_covid_pca ~ t_, data = may)
> m.covid.june <- lm(china_covid_pca ~ t_, data = june)
> m.covid.may.waves <- lm(china_covid_pca_waves ~ t_, data = may)
> m.covid.june.waves <- lm(china_covid_pca_waves ~ t_, data = june)
> 
> us.gov.may <- lm(usa_gov_pca ~ t_, data = may)
> us.gov.june <- lm(usa_gov_pca ~ t_, data = june)
> us.gov.may.waves <- lm(usa_gov_pca_waves ~ t_, data = may)
> us.gov.june.waves <- lm(usa_gov_pca_waves ~ t_, data = june)
> 
> us.people.may <- lm(usa_people_pca ~ t_, data = may)
> us.people.june <- lm(usa_people_pca ~ t_, data = june)
> us.people.may.waves <- lm(usa_people_pca_waves ~ t_, data = may)
> us.people.june.waves <- lm(usa_people_pca_waves ~ t_, data = june)
> 
> us.policy.may <- lm(usa_policy_pca ~ t_, data = may)
> us.policy.june <- lm(usa_policy_pca ~ t_, data = june)
> us.policy.may.waves <- lm(usa_policy_pca_waves ~ t_, data = may)
> us.policy.june.waves <- lm(usa_policy_pca_waves ~ t_, data = june)
> 
> us.covid.may <- lm(usa_covid_pca ~ t_, data = may)
> us.covid.june <- lm(usa_covid_pca ~ t_, data = june)
> us.covid.may.waves <- lm(usa_covid_pca_waves ~ t_, data = may)
> us.covid.june.waves <- lm(usa_covid_pca_waves ~ t_, data = june)
> 
> #Pooled data (not necessarily pooled PCA!)
> # PCA fit on pooled data
> m.gov <- lm(china_gov_pca ~ t_, data = pooled)
> # PCA fit separate for each wave
> m.gov.waves <- lm(china_gov_pca_waves ~t_, data = pooled)
> 
> m.people <- lm(china_people_pca ~ t_, data = pooled)
> m.people.waves <- lm(china_people_pca_waves ~t_, data = pooled)
> 
> m.policy <- lm(china_policy_pca ~ t_, data = pooled)
> m.policy.waves <- lm(china_policy_pca_waves ~t_, data = pooled)
> 
> m.covid <- lm(china_covid_pca ~ t_, data = pooled)
> m.covid.waves <- lm(china_covid_pca_waves ~t_, data = pooled)
> 
> #######################################################
> # Make figures
> ########################################################
> myoutcomes <- c("Government", "People", "Policy", "Covid-19")
> 
> #Fig 1 done with each loading scheme
> models <- list(m.gov, m.people, m.policy, m.covid,
>                m.gov.waves, m.people.waves, m.policy.waves, m.covid.waves)
> outcome <- myoutcomes
> term <- 2
> estimate <- unlist(lapply(models, function(x)coef(x)[[term]]))
> se <- unlist(lapply(models, function(x)coef(summary(x))[term, 2]))
> lower <- estimate - 1.96*se
> upper <- estimate + 1.96*se
> df <- cbind(outcome, estimate, se, lower, upper) %>% as.data.frame()
> df$outcome <- factor(df$outcome, levels = rev(outcome))
> df$estimate <- df$estimate %>% as.character() %>% as.numeric()
> df$lower <- df$lower %>% as.character() %>% as.numeric()
> df$upper <- df$upper %>% as.character() %>% as.numeric()
> df$loading <- c(rep("Pooled", 4), rep("By wave", 4))
> 
> fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome))
> fig1 <- fig1 + geom_point()
> fig1 <- fig1 + geom_errorbarh(height = 0)
> fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
> fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
>                                                           panel.grid.minor = element_blank(),
>                                                           panel.background = element_blank(),
>                                                           axis.title.y = element_blank(),
>                                                           legend.title = element_blank())
> fig1 <- fig1 + facet_wrap(~ loading) + theme(plot.title = element_text(hjust = 0.5))
> pdf("FigureA4.pdf", width = 4.375, height = 3.5)
> fig1
> dev.off()
> 
> # Fig 2 done with each loading scheme
> # In original, faceted by country.
> # This code will make a figure for each country, faceted by loading scheme.
> 
> # China
> term <- 3
> mymodels <- list(m.gov.may, m.gov.june,
>                  m.people.may, m.people.june,
>                  m.policy.may, m.policy.june,
>                  m.covid.may, m.covid.june,
>                  m.gov.may.waves, m.gov.june.waves,
>                  m.people.may.waves, m.people.june.waves,
>                  m.policy.may.waves, m.policy.june.waves,
>                  m.covid.may.waves, m.covid.june.waves)
> myoutcomes <- c("Government", "People", "Policy", "Covid-19")
> outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
> wave <- rep(c("normal", "crisis"), length(myoutcomes))
> estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
> se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
> lower <- estimate - 1.96*se
> upper <- estimate + 1.96*se
> df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
> df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
> df$wave <- factor(df$wave, levels = c("normal", "crisis"))
> df$estimate <- df$estimate %>% as.character() %>% as.numeric()
> df$lower <- df$lower %>% as.character() %>% as.numeric()
> df$upper <- df$upper %>% as.character() %>% as.numeric()
> df$loading <- c(rep("Pooled", 8), rep("By wave", 8))
> 
> fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
> fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
> fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
> fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
> fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
> fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
> fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
>                                                           panel.grid.minor = element_blank(),
>                                                           panel.background = element_blank(),
>                                                           axis.title.y = element_blank(),
>                                                           legend.title = element_blank())
> fig1 + facet_wrap(~ loading)
> 
> 
> pdf("FigureA5.pdf", width = 4.5, height = 3.5)
> fig1 + facet_wrap(~ loading)
> dev.off()
> 
> # USA
> term <- 3
> mymodels <- list(us.gov.may, us.gov.june,
>                  us.people.may, us.people.june,
>                  us.policy.may, us.policy.june,
>                  us.covid.may, us.covid.june,
>                  us.gov.may.waves, us.gov.june.waves,
>                  us.people.may.waves, us.people.june.waves,
>                  us.policy.may.waves, us.policy.june.waves,
>                  us.covid.may.waves, us.covid.june.waves)
> myoutcomes <- c("Government", "People", "Policy", "Covid-19")
> outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
> wave <- rep(c("normal", "crisis"), length(myoutcomes))
> estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
> se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
> lower <- estimate - 1.96*se
> upper <- estimate + 1.96*se
> df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
> df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
> df$wave <- factor(df$wave, levels = c("normal", "crisis"))
> df$estimate <- df$estimate %>% as.character() %>% as.numeric()
> df$lower <- df$lower %>% as.character() %>% as.numeric()
> df$upper <- df$upper %>% as.character() %>% as.numeric()
> df$loading <- c(rep("Pooled", 8), rep("By wave", 8))
> 
> fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
> fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
> fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
> fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
> fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
> fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
> fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
>                                                           panel.grid.minor = element_blank(),
>                                                           panel.background = element_blank(),
>                                                           axis.title.y = element_blank(),
>                                                           legend.title = element_blank())
> fig1 + facet_wrap(~ loading)
> 
> 
> pdf("FigureA6.pdf", width = 4.5, height = 3.5)
> fig1 + facet_wrap(~ loading)
> dev.off()
> 
> 
> # Fig 3 done with each loading scheme
> term <- 2
> mymodels <- list(m.gov.may, m.gov.june,
>                  m.people.may, m.people.june,
>                  m.policy.may, m.policy.june,
>                  m.covid.may, m.covid.june,
>                  m.gov.may.waves, m.gov.june.waves,
>                  m.people.may.waves, m.people.june.waves,
>                  m.policy.may.waves, m.policy.june.waves,
>                  m.covid.may.waves, m.covid.june.waves)
> myoutcomes <- c("Government", "People", "Policy", "Covid-19")
> outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
> wave <- rep(c("normal", "crisis"), length(myoutcomes))
> estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
> se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
> lower <- estimate - 1.96*se
> upper <- estimate + 1.96*se
> df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
> df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
> df$wave <- factor(df$wave, levels = c("normal", "crisis"))
> df$estimate <- df$estimate %>% as.character() %>% as.numeric()
> df$lower <- df$lower %>% as.character() %>% as.numeric()
> df$upper <- df$upper %>% as.character() %>% as.numeric()
> df$loading <- c(rep("Pooled", 8), rep("By wave", 8))
> 
> fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
> fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
> fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
> fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
> fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
> fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
> fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
>                                                           panel.grid.minor = element_blank(),
>                                                           panel.background = element_blank(),
>                                                           axis.title.y = element_blank(),
>                                                           legend.title = element_blank())
> fig1 + facet_wrap(~ loading)
> 
> 
> pdf("FigureA7.pdf", width = 4.5, height = 3.5)
> fig1 + facet_wrap(~ loading)
> dev.off()
> 
> log_close()

========================================================================= 
Log End Time: 2022-06-18 13:21:43 
Log Elapsed Time: 0 00:00:07 
========================================================================= 
